Scattering matrix of the boundary of a nonlocal metamaterial 
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We present a simple model of the interface between a local homogeneous medium and a potentially 
nonlocal metamaterial/photonic crystal. This model allows us to calculate the scattering matrix 
elements of the interface for a plane wave of light normally incident upon the interface from either 
direction. The resulting scattering matrix provides insight into the non-Maxwellian boundary condi- 
tions present at the interface between a homogeneous medium and a metamaterial/photonic crystal 
with strong spatial dispersion. We present the model mathematically. As an example, the model 
is used to calculate the scattering matrix of the interface between vacuum and a simple photonic 
crystal. Several tests of the calculated scattering matrix elements are presented. Finally, we used 
the results of the scattering model to postulate possible forms for the non-Maxwellian boundary 
conditions. 
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I. INTRODUCTION 

Most of the theoretical/numerical research into meta- 
materials since the inception of the fieldi has been fo- 
cused on designing and modelling new types of meta- 
material crystals and inclusions, as well as developing 
applications for metamaterials. A much smaller subset 
of research has investigated different methods of char- 
acterizing metamaterials using numerical simulations, a 
small number of notable examples being Refs. 0-i which 
characterize metamaterials (and metasurfaces) with scat- 
tering simulations as well as Refs. IB-fill which use non- 
scattering methods to attempt to homogenize metama- 
terials. Another small subset of metamaterial research, 
which is relevant to many of the metamaterial character- 
ization/homogenization methods, is the investigation of 
boundary conditions at the interface between metama- 
terials and other media. This research can be broadly 
divided into two categories. First, there are additional 
boundary conditions^—, which are extra boundary con- 
ditions in addition to the standard Maxwellian bound- 
ary conditions that are necessary at the interface of a 
medium because of the presence of multiple propagating 
modes due to spatial dispersion. Second, there are non- 
Maxwellian boundary conditions—'—, which are modifi- 
cations of the standard Maxwellian boundary conditions 
that become necessary in media with spatial dispersion, 
even in the absence of extra propagating modes. 

In this paper we only consider media supporting single 
propagating modes and as such are only concerned with 
non-Maxwellian boundary conditions. In particular, we 
ask the question: what is the scattering matrix for the 
interface between a homogeneous medium and a homog- 
enized metamaterial? To see why this is an interesting 
question, consider the example of the interface between 
two different homogeneous media and then contrast that 
with the metamaterial question. 

Assume we have two semi-infinite homogeneous me- 
dia, connected by a sharp interface, both described by 
isotropic, local (non-spatially dispersive), but potentially 
temporally dispersive constitutive parameters. The two 
media, which will be labelled medium 1 and medium 2, 



are described by the scalar permittivities e\ and 62, as 
well as the scalar permeabilities /zi and ^2, respectively. 
For plane waves normally incident upon the interface, 
the scattering at the interface depends on the r elative 
impedances of the two media, defined as z\ = [J-i/ti 
and Z2 = \J 1^2/^2- The scattering matrix elements them- 
selves will be labelled as for the electric field reflection 
amplitude of a plane wave incident upon the interface 
from medium i, and tij for the electric field transmission 
amplitude of a plane wave incident upon the interface 
from medium i and transmitted into medium j. Here 
and for the rest of the paper, all fields are assumed to be 
monochromatic. For our example of two homogeneous 
media, the formulas for the for scattering matrix elements 
are 



r%2 



T21 



*21 
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(1) 



Another convenient way to represent these scattering ma- 
trix elements, is to define the last three in terms of r\2 



T2\ = -ri2, 

t\2 = 1 + fl2, 

t 2 l = 1 + T21 = 1 - r i2 . 



(2) 



Eqs. (JTJ2J) are both results of the so called Maxwellian 
boundary conditions, which for normally incident waves 
is the requirement that the tangential electric (E) and 
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magnetic (H) fields be continuous across a boundary. For 
a metamaterial with strong spatial dispersion (or any ma- 
terial with strong spatial dispersion), we do not expect 
these boundary conditions to hold true. This is easy to 
see by considering the Poynting flux in a material with 
strong spatial dispersion. In such a material, a second 
term must be added the the standard Poynting flux-22r— . 
The total time averaged Poynting flux is 

S = S° + S 1 , (3) 

where S° is the standard expression for the time averaged 
Poynting flux 

s o = R^xtn (4) 

and S 1 is an extra term due to spatial dispersion 

si =-'H E, -s- E+E, -I- H 

(5) 

where is a unit vector in the i-th direction and there 
is an implicit sum over all values of i. Also, e and jl 
are the permittivity and permeability respectively and 
£ and C are the bianisotropic constitutive parameters^. 
Here and throughout the rest of this paper we are using 
Heaviside-Lorentz unitsSS, which are like Gaussian units 
with the 4-7T factor absorbed into the definition of the 
charges and currents. 

For a plane wave, normally incident upon an inter- 
face between a homogeneous medium and a metamaterial 
crystal, it is easy to see from Eqs (JS])-© that when strong 
spatial dispersion is present, the tangential electric and 
magnetic fields cannot both be continuous. This is be- 
cause continuous tangential electric and magnetic fields 
would force the component of S° normal to the interface 
to be continuous. If S 1 is non-zero due to the presence 
of spatial dispersion, this would prevent the total time 
averaged Poynting flux S from being continuous across 
the interface. Since the relations between the scattering 
matrix elements of the interface in Eq. @ are derived us- 
ing the Maxwellian boundary conditions, we can expect 
that in the presence of strong spatial dispersion some or 
all of the relations in Eq. ([2]) should fail. Understand- 
ing the boundary conditions of metamaterials with spa- 
tial dispersion is essential to a complete understanding 
of metamaterials. Calculating the scattering matrix of a 
metamaterial interface is an important first step towards 
understanding these new boundary conditions. The pur- 
pose of this paper is to present a simple model for nu- 
merically calculating these scattering matrix elements. 



The authors are only aware of one previous attempt to 
calculate the scattering matrix of a metamaterial inter- 
face^. This method involves inferring the interface scat- 
tering matrix elements as well as the index of refraction 
of a metamaterial from the total reflection and transmis- 
sion amplitudes for two different slabs of the same meta- 
material with different thicknesses. While the authors of 
Ref. dl| acknowledge that the Maxwellian boundary con- 
ditions can fail at a metamaterial boundary, they make 
the assumption that the two transmission coefficients of 
the metamaterial interface, labelled i 12 and t 2 \ in this 
paper, are equal due to Lorentz reciprocity. This leaves 
only three scattering matrix elements plus the index of 
refraction to be algebraically inferred from the four total 
scattering matrix elements of the two slabs of differing 
thickness. We disagree with this assumption and point 
out that it is not true even for the interface between vac- 
uum and a simple dielectric, as can be seen from Eq. ([1}. 

We also note that it is well known that one can use Rig- 
orous Coupled Wave Analysis^ (also known as Fourier 
Modal Method) to calculate the scattering matrix be- 
tween plane waves in a vacuum (or some other homo- 
geneous local medium) and the Bloch modes of a meta- 
material/photonic crystal, separated by a sharp inter- 
face. While this method can correctly calculate the re- 
flection coefficient back into vacuum, labelled r\i in this 
paper, the remaining scattering amplitudes (plane wave 
to Bloch mode, Bloch mode to plane wave, and Bloch 
mode to Bloch mode) must not be confused with the 
scattering matrix elements calculated in this paper. The 
scattering matrix elements in this paper are the ratios 
of electric field amplitudes of planes waves, both in the 
homogeneous medium and in the homogenized metama- 
terial/photonic crystal. 

In Sec. |TT] we present the numerical model for calcu- 
lating the scattering matrix of the metamaterial inter- 
face. This model is based on a similar model used for 
metamaterial homogenization^. In addition, this scatter- 
ing model relies upon constitutive parameters calculated 
with the same homogenization model. In Sec. IIIII we ap- 
ply this scattering model to calculate the scattering ma- 
trix of the interface between vacuum and a simple pho- 
tonic crystal. We then present two tests of the calculated 
scattering matrix elements. Finally, in Sec. IIV1 we use 
the results of the interface scattering model to postulate 
a set of phenomenological boundary conditions. 



II. SCATTERING MATRIX OF THE 
INTERFACE BETWEEN A HOMOGENEOUS 
MEDIUM AND A METAMATERIAL 

A. The constitutive parameters of a one 
dimensional metamaterial 

The one dimensional model that we present in this 
paper for calculating the scattering matrix of a meta- 
material interface is strongly related to the one dimen- 
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sional model used for metamaterial homogenization in 
Ref. ID. The constitutive parameters used in this paper 
are calculated using the same homogenization method, 
originally presented for p-polarized waves in Ref. 1, and 
re-presented for s-polarized waves (E = E z z and H = 
H^x + H y y) in the appendix of this paper. Here we 
briefly review the properties of the constitutive parame- 
ters determined from the one dimensional metamaterial 
homogenization model for s-polarized waves. 

In this paper we will only consider two dimen- 
sional highly symmetric metamaterial/photonic crystals, 
though the concept is easily generalized to highly sym- 
metric three dimensional crystals. In two and three di- 
mensions the symmetry requirements are that the crys- 
tal have a rectangular unit cell, and that the crystal has 
reflection symmetry in two directions tangential to the 
interface, making the direction normal to the interface a 
principle axis of the crystal. In three dimensions, there 
is also the possibility of homogenizing a crystal with a 
rectangular unit cell as well as 180° rotational symme- 
try around the direction normal to the interface, again 
making this direction a principle axis, even without re- 
flection symmetry transverse to the interface. In two 
and three dimensions, reflection symmetry in the direc- 
tion of propagation (normal to the interface) is not neces- 
sary. Also, for our limited two dimensional example, we 
are calculating scattering matrix elements for s-polarized 
waves, though it is straightforward to do the same for 
p-polarized waves. 

The constitutive relationship for a s-polarized electro- 
magnetic field that is harmonic in space and time, prop- 
agating in the x-direction in a two dimensional crystal 
is 



\ B y) \ H y/ 
where C is the constitutive matrix 



(6) 



C(u,k x ) = 



(e zz {w,k x ) £, zy {uj,k x )\ 



(7) 



We emphasize that all of the constitutive parameters are 
functions of frequency u and wavenumber k x , the de- 
pendence on the wavenumber determining the degree of 
spatial dispersion, which in turn is the cause of the non- 
Maxwellian boundary conditions. 

As shown in Ref. if the metamaterial/photonic crys- 
tal of interest is reciprocal, the constitutive parameters 
obey the Lorentz reciprocity relations for a spatially dis- 
persive material 



e zz (cj,-k x ) = e zz (u,k x ), 

Hyy(u},-kx) = fJ,yy(u),k x ), 
£ Z y(u,-k x ) = -Cy Z (u,k x ). 



(8) 



This allows us to represent the constitutive parameters 
as 



C = 



[lyy 



(9) 



Here e yy , (j, zz and K e are even functions of k x but k q is an 
odd function of k x . In this paper we will only consider 
reciprocal metamaterial crystals. 

From the constitutive parameters we can define an 
impedance for freely propagating waves 



z ± = T 



k x (oj) + L0K o / C =F UK e /c 
k x (uj) + u>Ko/c ± LOK e /c 



(10) 



ue zz /c 



Here z is the impedance for a plane wave propagating 
in the positive (+) or negative (-) x-direction, the dif- 
ference between the two due to possible asymmetry of 
the crystal in the direction of propagation and the re- 
sulting nonzero n e . The wavenumber k x (uj) in Eq. (|10[) 
is the wavenumber of a freely propagating eigenmode of 
the crystal propagating in the positive x-direction. For a 
passive crystal this implies that Im(k x ) < assuming the 
convention that harmonic plane waves vary as ^{ut-kxx) _ 
The constitutive parameters in Eq. ([TO)) are also evalu- 
ated using this positively propagating k x (ui). 

Finally, for a crystal with symmetry in the direction of 
propagation (here the x direction), symmetry forces n e — 
0. Spatial dispersion however, still allows for nonzero 
bianisotropy through k q . 



B. Scattering model of the metamaterial interface 

The model presented in this section for calculating the 
scattering matrix of a metamaterial interface builds on 
a earlier model briefly described in Ref. HI for calculat- 
ing the reflection from a semi-infinite metamaterial. This 
newer, more complete version of the model allows for the 
calculation of the entire scattering matrix of the metama- 
terial interface. Fig. Q] provides a diagram of the meta- 
material interface. The incident waves are normally in- 
cident upon the interface, the normal direction being de- 
fined as the x-direction. As before, the scattering matrix 
elements are labelled as r,j for electric field reflection am- 
plitude of a plane wave incident upon the interface from 
medium i, and tij for the electric field transmission am- 
plitude of a plane wave incident upon the interface from 
medium i and transmitted into medium j. The homoge- 
neous medium is labelled medium 1 and the homogenized 
metamaterial is labelled medium 2. 

Just as in Ref. [^, we model the metamaterial crystal 
by replacing individual layers of the crystal with equiva- 
lent metasurfaces characterized by a surface polarizabil- 
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Homogeneous M p tamat p ria | 
Medium Metamatenai 




FIG. 1: A diagram of the scattering problem examined in 
this paper. The scattering occurs at the interface between a 
semi-infinite homogeneous medium and a semi-infinte meta- 
material. The metamaterial is potentially nonlocal (spatially 
dispersive). 



ity. The metasurfaces only interact with adjacent meta- 
surfaces through planes waves. Since all evanescent inter- 
actions are ignored, higher order propagating modes are 
absent from the model and the corresponding additional 
boundary conditions are not supported in the model. As 
m Ref.i, the one dimensional model of the interface has a 
algebraic solution which can easily be solved numerically. 
There are six degrees of freedom in the model. These 
are the coefficients and bi, shown in Fig. [3J which are 
the electric field amplitudes for plane waves travelling 
in the positive x-direction (aj) and negative x-direction 
(pi). The microscopic electric and magnetic fields at any 
point in the domain of our one dimensional model are 



-ikiX 



-ikiX 



(11) 
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Zb 



(12) 

Next, we consider the relationship between the field co- 
efficients of the second and third subdomains. These two 
subdomains are separated by a metasurface. This is an 
infinitesimally thin layer with the boundary condition 



(13) 




Here, Ah,, 



are the mi- 



±y = h+ — h~, where h+ and h~ 
croscopic magnetic fields on the right and left sides of 
the metasurface respectively. Similarly, Ac z = e+ — c~ 
where ej and e~ are the microscopic electric fields on 
the right and left sides of the metasurface respectively. 
The local fields are defined as E' 



lor 



8+ + ej)/2 and 

(h+ + hy)/2. Finally, the electromagnetic re- 
sponse of the metasurface is characterized by the surface 
polarizability 



a(ui) 



'zy 



(14) 



//</ 



which, can be derived from the numerically calculated 
scattering matrix of a single layer of the metasurface^ 



Here and bi are the appropriate field coefficients for 
the particular value of x. Also, = v/^Ai is the 
impedance and ki 



Also, Zi 

/€ifiiU>/c is the wavenumber for 
the corresponding subdomain. In general, the constitu- 
tive parameters of the homogeneous medium t\ and \i\ 
can be different from the constitutive parameters of the 
background (substrate) of the metamaterial e,; and /z, for 
i = 2,3. 

Only a single metasurface is necessary to represent 
the semi-infinite array. It is trivial to include additional 
metasurfaces representing additional layers of the crys- 
tal, but this is unnecessary since it results in the same 
final interface scattering matrix. Since there are six de- 
grees of freedom, we require six equations of constraint. 
These include four boundary conditions at the two inte- 
rior boundaries (two for each boundary), as well as one 
boundary condition at each of the two exterior bound- 
aries. The first interior boundary to consider is the one 
between the homogeneous medium and the background 
medium (substrate) of the metamaterial located at x = 0. 
Due to Maxwellian boundary conditions at the interface 
between the vacuum and the background material (conti- 
nuity of tangential electric and magnetic fields), the rela- 
tionship between the 1-st and 2-nd sets of field coefficients 
is 




FIG. 2: A diagram of the simple one dimensional model used 
to solve the metamaterial interface scattering problem. There 
are three domains, each with two s-polarized plane waves, one 
propagating in the positive x direction with electric field am- 
plitude ai, and one propagating in the negative x direction 
with electric field amplitude The 1-st and 2-nd domains are 
separated by an interface between the homogeneous medium 
and the background (substrate) medium of the metamaterial. 
The 2-nd and 3-rd domains are separated by a metasurface 
which serves as a subsitute for a single layer of the metama- 
terial crystal one unit cell thick. 
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f[l + det(S)-(S 11 + S 22 )]/z b (Si 



S 2 



(5*12 — S21) 



u/c(l + S 12 + S21 - det(S)) I -(S n - S 22 ) - (S12 ~ S 21 ) [1 + dct(S) + {S n + S 22 )} 



(15) 



The scattering matrix of the metasurface in Eq. (fT5)l is defined with the reference plane centered at the location of 
the metasurface. We also note that Ref. ID presents a similar equation for the surface polarizability of a metasurface 
interacting with a p-polarized wave, as well as instructions for how to transform the p-polarizcd equation into an 
s-polarized version. These, instructions were incorrect. Eq. (| 15|) is the correct equation for the surface polarizability 
of a metasurface interacting with a s-polarized wave. 

Using Eq. (fl3|) . we can relate the difference in the fields across the metasurface to the field coefficients with 



Mh, 
lAe, 



= z b 



3 — ifei,a/2 



ikbCi/2 
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ik b a/2 
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Using the definition of the local fields, we can relate the local fields to the field coefficients with 



TjloC 



( £ -ifc k o/2 



± e ik b a/2 fg-ifci.o/2 



-ik^a/2 
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2z~b 
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2z b 
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2V 



ik b a/2 
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J 


W 



M 2 



(16) 



(17) 



Thus, the two boundary conditions in Eq. (1131) can be represented as 



Mi - i—a ■ M 2 
c 



«3 



(18) 



If the fields associated with the field coefficients in Eq. (fT5)) are forced to be Bloch periodic, that is if the fields at 
x = a differ from the fields at x — by the Bloch amplitude e _lfex<1 , then the Bloch periodicity combined with the 
boundary conditions in Eq. (|18j) allow us to derive a dispersion relation for the one dimensional array of metasurfaces 



2 a ee mm _ era me s. me , era 

c z 4 / c 2 



(19) 

I - -- —yy - ~*v-y* I cos(k a) - ^ sm(fc a). 



This dispersion relation for s-polarized free waves prop- 
agating in a one dimensional array of metasurfaces is 
electromagnetically dual to the dispersion relation for p- 
polarized waves presented in Ref. |9|. 

We have now imposed four of the six required equa- 
tions of constraint. There remain two exterior boundary 
conditions that must be imposed in order to solve for the 
six field coefficients. These two exterior boundary condi- 



tions involve the electric field amplitudes of plane waves 
that are incident upon the metamaterial interface. First 
we consider a plane wave incident from vacuum. The 
constraint for this is simply 

a x = E^, (20) 

where E™ c is the electric field amplitude of the plane 
wave incident upon the interface from the homogeneous 
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medium. The constraint for a plane wave incident from 
the metamaterial is more complicated. Eq. (IT71) relates 
the field coefficients to the local fields, and the relation 9 



Y- 




I or 



Y 



relates the local fields at the metasurface to E! 



(21) 
, the 



electric field amplitude of the plane wave incident upon 
the interface from the metamaterial and E^ 1 *, the elec- 
tric field amplitude of the plane wave propagating in the 
metamaterial outwards from the interface. Both E™ c 
and EJ^' are electric field amplitudes at the location of 
the interface, x is the macroscopic susceptibility de- 
fined in Eq. (|A13I) for plane waves moving in the pos- 
itive (x+ = x(u,k x )) and negative (x - = x(w, -fc x )) 
x-directions, where k x {u>) is the wavenumber of a freely 
propagating eigenmode of the crystal moving in the pos- 
itive x-direction determined from Eq. (|19|) . Also, z ± is 
the impedance of the homogenized metamaterial defined 
in Eq. (I10[) calculated with wavenumber k x (oj) provided 
by the dispersion relation Eq. (JTHJ). Both the left and 
right hand sides of Eq. (|2ip are equal to the polarization 
density of the homogenized metamaterial at the location 
of the metasurface. 

Combining Eqs. (fTTj) and ([2~Tj) . we get 
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(22) 



This provides us with our final equation of constraint 



This small algebraic system of equations is easily solved 
numerically. 

In order to calculate the scattering matrix for the 
metamaterial interface, we must solve Eq. (I24[) twice, 
once with a plane wave incident from the homogeneous 
medium, and once with a plane wave incident from 
the metamaterial. First we solve the case of a plane 
wave incident from the homogeneous medium by setting 
E™ c = and setting E™ c to an arbitrary non-zero value. 
After solving Eq. (j2"4")l for the field coefficients, we can 
calculate two of the scattering matrix elements of the 
interface, 



r\2 



E'. 



ho = 



El 



:(!,())■ X- 



b 2 



(25) 



Next, we solve the case of a plane wave incident from 
the metamaterial by setting E™ c = and setting E™ c to 
an arbitrary non-zero value. After solving Eq. (|24|) for 
the field coefficients, the remaining two scattering matrix 
elements are 



1 



r 2 l 



hi 



•pine 

h 



(1,0) -X- 



b 2 



(26) 



III. 



EXAMPLE: SCATTERING MATRIX OF A 
PHOTONIC CRYSTAL INTERFACE 



(0,1) -X- 



(a 2 \ 
b 2 



E! 



(23) 



We now have six equations of constraint to allow us to 
solve our model. Putting them all (Eqs. (fHj) . (IT51) . ([2U|) 
and (|23p ) together, they are 
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(24) 



A. Dispersion relation and scattering matrix 

As a demonstration of the metamaterial interface scat- 
tering model, we calculate the scattering matrix of the 
interface between vacuum and a simple two dimensional 
photonic crystal, the unit cell of which is shown in Fig. [3] 
The unit cell is square with lattice constant a. At the cen- 
ter of the unit cell is a cylinder with radius r — 0.3a. The 
cylinder is a dielectric with permittivity e = 12 — i • 10~ 3 
and is surrounded by vacuum with permittivity e = 1. 
The negative imaginary part in the permittivity of the 
cylinder indicates loss. Fig.[3]also shows a complex k x (ui) 
dispersion relation for a one dimensional array of meta- 
surfaces corresponding to the cylindrical photonic crystal 
and calculated with Eq. (|T9l) . Notice there are two band 
gaps, indicated by the non-zero imaginary part of k x and 
shaded in Fig. [3] and in the remaining figures. This dis- 
persion relation agrees very well with the complex k x (uj) 
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0.5 1 1.5 2 2.5 3 



wa/c 

FIG. 3: (a) Diagram of the unit cell of the photonic crystal ex- 
amined in this section, (b) Complex k x dispersion relation for 
an s-polarized eigenmode propagating along a principle axis 
of the photonic crystal. Band gaps are indicated by nonzero 
\m(k x ) and are shaded. 



dispersion relation of the cylindrical photonic crystal cal- 
culated from a finite element simulation— (not shown), 
except for a narrow band near uj = 2.5c/ a. 

The constitutive parameters of the photonic crystal 
must first be determined if we are to calculate the scat- 
tering matrix of the interface. This is done with the 
procedure presented for p-polarized waves in Ref . i, and 
re-presented in the appendix for s-polarized waves. First, 
the surface polarizability of a single layer of the crystal 
is calculated from the scattering matrix of a single layer. 
This is then used to calculate a set of macroscopic con- 
stitutive parameters relating the electric and magnetic 
fields to the electric displacement and magnetic flux den- 
sity fields 

( D 'W £ " ,27, 

\ByJ \K llyyj \RyJ 

Each of these constitutive parameters is a function of 
both lu and k x , the dependence on k x indicating spa- 
tial dispersion. Since we are only considering s-polarized 
waves, and due to the symmetry and reciprocity of the 
crystal unit cell, we only need to consider three consti- 
tutive parameters. Also, as mentioned in Sec. flTR. be- 
cause of the reciprocity of the crystal, both e zz and \i yy 
are even functions of k x , and k is an odd function of 
k x . Using these constitutive parameters, along with the 
model described in the previous section, we can calculate 
the scattering matrix for the interface between vacuum 




0.5 1 1.5 2 2.5 3 

toa/c 



FIG. 4: (a) Absolute value and (b) argument of the reflection 
amplitude rfj™ calculated from a finite element simulation as 
well as the scattering matrix elements m, r^\, t\2 and £21 of 
the vacuum/crystal interface calculated from the model pre- 
sented in Sec.HTB. (c) Test of the relations in Eq. (f2]). Shaded 
regions indicate a band gap. The relation ri2 = — J"2i appears 
to be valid while the relations in Eq. ((2} for the transmission 
coefficients fail. Notice that the scattering matrix elements 
are real valued outside the band gap but complex valued in- 
side the band gap. Also, in the band gap | ri2 1 = [fai| = 1. 

and the photonic crystal. The four parameters of the 
scattering matrix are plotted in Fig. SJ 

The first important result in Fig. [4] is the agreement 
between r{| m , which is the the reflection amplitude of a 
semi-infinite photonic crystal calculated from a finite ele- 
ment simulation and r\i , which is a scattering matrix el- 
ement calculated from the model presented in Sec-HIB. It 
is not known how to calculate the other scattering matrix 
elements from a finite element simulation and therefore 
it is not possible to test them in this way. The second 
result, is that the two reflection coefficients r\% and r 2 i 
do in fact appear to be the negative of each other. This 
can be seen more clearly in Fig. [4}:, which plots the values 
\r%\ +ri2|, |ti2 — 1 — rial and |t 2 i - 1 + ri2|. If the rela- 
tions in Eq. ^ are correct then the expressions plotted 
in Fig. 0J: should all be zero. We see from Fig. 0h that 
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FIG. 5: Comparision of reflection and transmission amplitudes for a photonic crystal slab three units cell thick, calculated 
in three ways. From a finite element simulation (R fem and T /em ), from Eq. ([28]) using all four scattering matrix elements 
calculated from the model (R and T), and from the simplified expressions in Eq. Q29p calculated from ri2 (Rr 12 an d T ri2 ). 
(a) Absolute value and (b) argument of the reflection amplitudes, (c) Absolute value and (d) argument of the transmission 
amplitudes. Shaded regions indicate a band gap. Note that the apparent n jumps in arg(R) occur when |R| approaches zero. 



this indeed true for the relation ri\ = —T\i. However, 
the two remaining relations in Eq. ([2]), while valid in the 
long wavelength limit, do appear to fail as expected due 
to spatial dispersion. 



not used any of the relations in Eq. ^ which are often 
used to simplify these types of expressions. If we had 
used Eq. ©, the reflection and transmission amplitudes 
could have been expressed as 



B. Test: Transmission and reflection of a photonic 
crystal slab 

As a test of the calculated interface scattering matrix, 
we can use the scattering matrix elements to calculate 
the reflection and transmission through a photonic crys- 
tal slab. The calculated reflection and transmission am- 
plitudes can then be compared against values calculated 
from a finite element simulation of the photonic crys- 
tal slab. For a photonic crystal slab, three layers thick 
and with two boundaries with vacuum, each boundary 
described by scattering matrix elements calculated from 
our model, the total reflection and transmission ampli- 
tudes for the slab are 



-6iA.-. r 



R 



ri2 + 



*i2r2i*2ie 



-&ik x a 



1 



£i2*2ie 



(28) 



1 



-Q\k x a ' 



Here k x (u) is a complex valued wavenumber calculated 
from Eq. (fH))) and plotted in Fig. [3] Notice that we have 



R ri2 
T, 12 



ri2 



1 ». 2 „— Qik x a 1 

i r 12 t 



(29) 



-3ik x a 



12 c 



In Fig. [5] we have plotted the complex valued reflec- 
tion amplitudes R^ em , R and R ri2 , and the transmission 
amplitudes T^ em , T and T ri2 . Here R fem and T^ em are 
calculated from a finite element simulation of a three lay- 
ered photonic crystal slab. 

We see from Fig. [5] that the transmission and reflection 
amplitudes calculated using Eq. (|28|) with the scattering 
matrix elements of the vacuum-crystal interface calcu- 
lated from our model agree very well with the finite ele- 
ment simulation. However, the simplified reflection and 
transmission amplitudes of Eq. (|29p using only ri2 also 
agree with the finite element simulation. This is a sur- 
prise since we have seen in Fig. 3] that with the exception 
of r2i = — ri2, the remaining relations of Eq. ([2j fail due 
to spatial dispersion. Comparing Eqs. (|28H29[) , we find a 
new relation between the scattering matrix elements, 



fi2*2i = (1 + ri 2 )(l + r 2 i), 



(30) 
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which is easily confirmed by directly examining the scat- 
tering matrix elements returned by our model. This new 
relation is true despite the fact that the relations for the 
transmission amplitudes in Eq. ^ fail due to spatial dis- 
persion. 

C. Test: Poynting flux at the vacuum-crystal 
interface 

There is an additional test we can perform to test the 
accuracy of the calculated scattering matrix elements of 
the vacuum-crystal interface. Energy flux at the interface 
must be conserved. If we use the calculated scattering 
matrix elements to calculate the macroscopic energy flux 
and compare it to the true microscopic energy flux, this 
provides us with a strong test of the scattering matrix 
elements. At the same time, calculating the energy flux 
is a good test of the nonlocal homogenization method 
we have used 9 to calculate the constitutive parameters 
of the photonic crystal. This homogenization procedure 
deliberately calculates constitutive parameters that are 
functions of the wavenumber. This is essential for cor- 
rectly calculating the Poynting flux since the term S 1 in 
Eq. ([5]) involves partial derivatives of the constitutive pa- 
rameters with respect to the wavenumber. To the best 
knowledge of the authors, there has been only one other 
peer reviewed paper where the nonlocal Poynting flux has 
been calculated inside of a metamaterial 30 . 

We can calculate the Poynting flux at the vacuum- 
crystal interface for two different scenarios. First, where 
a plane wave is incident upon the interface from the vac- 
uum and second, when a plane wave is incident upon the 
interface from the crystal. In the first case where the 
plane wave is incident from the vacuum, the macroscopic 
Poynting flux in the x direction is 



crystal, one incident and one reflected. For the case of 
two overlapping plane waves, each wave has a correspond- 
ing Si (Eq. ([5])) correction to the Poynting flux. This re- 
quires that the derivatives of the constitutive parameters 
be evaluated for both positive k x , for the reflected wave 
moving in the positive x-direction, as well as evaluated 
for negative k x , for the incident wave moving in the nega- 
tive x-direction. Assuming Lorentz reciprocity, ensuring 
that e zz , fj, yy and n e are even in k x and k q is odd in k x , 
we can relate the derivatives of the constitutive parame- 
ters evaluated at negative k x , to the derivatives evaluated 
at positive k x . Done correctly, the Poynting flux in the 
x direction at the vacuum-crystal interface with a plane 
wave incident from the crystal is 




It should be noted that the derivation for Eq. ([5p2r— 
assumes that the frequency u) and wavenumber k x are 
approximately real valued. This is not necessarily true 
for an evanescent wave or for a propagating wave in the 
presence of high losses. In either case the expression for 
the Poynting flux in Eqs. pij|33[) is in general not correct. 

The two Poynting fluxes S„ and S m are plotted in Fig|6l 
along with the microscopic Poynting flux in the x direc- 
tion for each case 



$v — — I— S^, , 

Re[E z H*] 



where 



= -Re 
2 



-Re [3+] M 2 K" 



77± — 



de 
dk x 

=p2ilm 



2Re 



\h2\ 2 \K n 



dn 
dk x 



dn e 1 9/x 
dk x \z ± \ 2 dk x 



(31) 



(32) 



Note that the derivatives of the constitutive parameters 
are evaluated for real valued wavenumbers k x . 

In the second scenario, with a plane wave incident from 
the metamaterial crystal, there are now two waves in the 



;Re[e z ig. 



(34) 



By definition, the microscopic fields at the interface agree 
with the macroscopic fields on the vacuum side of the 
interface. In both cases, the incident electric field is nor- 
malized to unity. We can see from Fig.|6l that the Poynt- 
ing flux calculated from the scattering matrix elements 
of the metamaterial interface and the constitutive pa- 
rameters of the homogenized photonic crystal agrees very 
well with the microscopically calculated Poynting flux for 
frequencies that lie in the passband. For frequencies in 
the bandgap, there is very little agreement between the 
macroscopic and microscopic Poynting fluxes. As men- 
tioned earlier, this is due to the large imaginary part of 
wavenumber k x of an evanescent wave in a band gap. 
However, the fact that we get the correct Poynting flux 
in the passband is positive evidence that the calculated 
scattering matrix elements are correct in the pass band. 

Both tests of the calculated interface scattering ma- 
trix elements returned positive results, however it should 
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FIG. 6: Poynting fluxes in a spatialy dispersive medium (ho- 
mogenized photonic crystal) at the metamaterial interface for 
(a) plane wave incident from the vacuum and (b) plane wave 
incident from the photonic crystal. Shaded regions indicate a 
band gap. In both cases the incident electric field is normal- 
ized to unity. We see good agreement between the spatially 
dispersive Poynting flux expressions S = S° + S 1 and the mi- 
croscopic Poynting flux S mlc for frequencies outside the band 
gap. 



be noted that the first test in Sec. IIIIB really indi- 
cates the validity of the relationships r%x = —r\2 and 
ii2*2i = (1 + r i2)(l + J"2i) between the various scat- 
tering matrix elements. Similarly, the test presented in 
Sec. IIIIC only tests the absolute value of the scattering 
matrix elements, and only provides positive results in the 
passbands since Eqs. (|3HI33|) are only valid in the pass- 
bands. 



IV. PHENOMENOLOGICAL 
NON-MAXWELLIAN BOUNDARY CONDITIONS 

A. Interface between local homogeneous medium 
and a symmetric crystal 

A surprising result of our investigation of the scattering 
matrix of the photonic crystal interface in Fig. [3]was that 
the relation ri2 = —^21 remains true even when spatial 
dispersion forces the remaining relations of Eq. ([2} to 
fail. Our test of the total reflection and transmission 
amplitudes for a slab of the cylindrical metamaterial led 
to the discovery of another relation in Eq. (|50"|) . 

In an attempt to try to generalize these relations, 
and use them to infer some insights into non-Maxwellian 
boundary conditions for metamaterials with spatial dis- 



Homogeneous Metamaterial I Homogeneous 



ri2 



*12 r 23 
r 2 i 



•r 3 2 



FIG. 7: Diagram of a slab of a nonlocal homogenized meta- 
material surrounded on both sides with semi-infinite local ho- 
mogeneous media with differing constitutive parameters. 



persion, we consider measuring the total reflection and 
transmission through a metamaterial slab, only this time 
the homogeneous media on the opposite sides of the slab 
have different constitutive parameters than vacuum and 
from each other. This simple geometry is shown in Fig. [7] 
Assuming that the three different media are lossless 
and that the metamaterial has reflection symmetry in the 
x-direction (n e — 0), and by performing the simple an- 
alytic derivation of the total reflection and transmission 
amplitudes of the slab, it is easy to see that conservation 
of energy can be guaranteed if the following relations are 
true 



1 = 1 ( l ~ r ™ 

z\ \1 + ri2 J 23 V 1 + r 32 
£12^23 = (1 + ri 2 ) (1 + r 2 z) 



(35) 



These relations can be easily confirmed numerically with 
the model presented in Sec. HTt3. 

The relation r\% = —r 2 i, along with Eqs. (|30I35[) allow 
us to guess the form of the boundary conditions for the 
interface between a local homogeneous medium and a 
nonlocal homogenized metamaterial. For the interface in 
Fig. [7] between media 1 and 2, the boundary conditions 
are 



1 + 7-12 = 




Ei 
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(36) 
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Here E*, E^, and are the electric and magnetic 
field amplitudes in media 1 and 2, and in every case 
the incident electric field is normalized to 1. The co- 
efficient a and b are interface parameters characteriz- 
ing the non-Maxwellian boundary condition at the in- 
terface. These phenomenological boundary conditions 
are also easily confirmed numerically with the scattering 
model presented in Sec. HlB. both for lossless and lossy 
metamaterials. The interface parameters a and b for the 
cylindrical photonic crystal shown in Fig. [3] are plotted 
in Fig. H 
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FIG. 8: The interface parameters of Eq. (|36|) for the nonlo- 
cal medium corresponding to the cylindrical photonic crystal 
shown in Fig. [3] 

There are a number of conclusions we can draw from 
Eq. (1551) . First, as can be seen in Fig. [FJ in the long 
wavelength limit, both a and b equal 1, corresponding 
to the standard Maxwellian boundary conditions. Sec- 
ond, the discontinuity in the transverse electric (E 2 ) and 
magnetic (H y ) fields is always proportional to the elec- 
tromagnetic field strength on the metamaterial side of 
the interface. The field strength on the side of the inter- 
face with the local homogeneous medium has no effect 
on the field discontinuities. Furthermore, the disconti- 
nuity in the magnetic field, which is usually associated 
with an electric surface current, is proportional to the 
magnetic field. Likewise, the discontinuity in the electric 
field, normally associated with a magnetic surface cur- 
rent, is proportional to the electric field. If one tried to 
characterize the metamaterial interface as a metasurface, 
as was attempted in Ref. [2l|, according to Eq. (j3"o) . the 
surface currents would be proportional to the electromag- 
netic fields on the metamaterial side of the interface, and 
the surface polarizability matrix would have zeros along 
its diagonal and nonzero values for its off diagonal com- 



ponents. Lastly, though this is not evident from Eq. (|36l) . 
it can be demonstrated numerically using the model pre- 
sented in Sec. UH 3 that the interface parameters a and 
6 only depend on the particular metamaterial. They are 
unchanged when the medium bordering the metamaterial 
is changed. We will return to this point in Sec. II VB . 



B. Interface between local homogeneous medium 
and an asymmetric crystal 

So far we have only considered crystals that have re- 
flection symmetry in the direction of wave propagation. 
We now briefly consider a crystal that does not have re- 
flection symmetry in the propagation direction, though 
it does still have reflection symmetries in the directions 
perpendicular to propagation. We cannot expect the re- 
lations we postulated in Eq. (|36p for a symmetric crystal 
to be correct for an asymmetric crystal. It is easy to see 
that in the long wavelength limit, we cannot even expect 
the relation ri2 ~ — 7"2i to be valid for an asymmetric 
crystal, due to the fact that in an asymmetric crystal the 
impedance depends on weather the electromagnetic wave 
is propagating in the positive or negative x-direction (see 
Eq. (TO])). 

There are a few thing we can assume about the bound- 
ary conditions for an asymmetric crystal. First, in the 
limit of an asymmetric crystal becoming symmetric, n e 
vanishes, z + — z~ 1 and the boundary conditions should 
reduce to Eq. ([36l) . Second, by examining Eqs. (|31ti33[) . 
we can see that any coefficients determining the bound- 
ary conditions for the interface with an asymmetric crys- 
tal, analogous to a and b for a symmetric crystal, should 
depend on the direction of propagation of the electro- 
magnetic waves. With these insights, we can guess the 
form of the non-Maxwellian boundary conditions for the 
interface of an asymmetric crystal. For an asymmetric 
crystal, that when homogenized becomes nonlocal and is 
labelled as medium 2, surrounded by two local homoge- 
neous media labelled 1 and 3, the boundary conditions 
for the interface between media 1 and 2 are 




a + r 2 i + a 

a +E 2 _+ a-E 2 z - 



(37) 
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and the boundary conditions for the interface between 
media 2 and 3 are 




(38) 



Here Z\ and Z3 are the impedances of media 1 and 3, and 
z 2 are the impedances of the asymmetric homogenized 
metamaterial given by Eq. (JTDJ). E^, E^ and E^ and 
H**, H2± and H3± are the electric and magnetic field 
amplitudes in media 1,2 and 3 for waves propagating in 
a positive (+) or negative (— ) x-directions, and in every 
case the incident electric field is normalized to 1 . Finally, 
and are interface parameters for waves propagating 
in the positive (+) and negative (— ) x-directions. 

As with the interface parameters for the symmetric 
metamaterial, the parameters and depend solely 
on the particular metamaterial and are not affected by 
the medium on the opposite side of the interface. Simi- 
larly, the discontinuity in the tangential electric and mag- 
netic fields is entirely due to the electromagnetic field 
strength on the metamaterial side of the interface. This 
can all easily be confirmed from the interface scattering 
model presented in Sec. HTt3. 

Finally, though the interface scattering model pre- 
sented in Sec. HlB is for an interface between a local 
homogeneous medium and a nonlocal homogenized meta- 
material, this model can easily be modified to calculate 
the scattering matrix elements of the interface between 
two different nonlocal homogenized metamaterials. The 
scattering matrix elements behave according to boundary 
conditions similar to Eqs. (|36H38[) . In this case the discon- 
tinuities in the tangential electric and magnetic fields are 
proportional to the electromagnetic fields on both sides of 
the interface. However, the boundary conditions are still 
determined by the interface parameters in Eqs. (|36H38p . 
That is to say, if one calculates the interface parameters 
for two different homogenized metamaterials using the 
scattering model presented in Sec. HlB. those scattering 
matrix elements can predict the scattering of an inter- 
face between those two different nonlocal homogenized 
metamaterials. 



V. CONCLUSION 

We have presented a simple model for calculating the 
scattering matrix for waves scattered from the inter- 
face between a homogeneous medium and a homogenized 
highly symmetric metamaterial crystal. In doing so we 
have verified that the relations in Eq. ([2]) between the var- 
ious scattering matrix elements are in general not true 
when spatial dispersion is present. Using this scatter- 
ing model, we have postulated a set of phcnomcnologi- 
cal boundary conditions for the interface of a metamate- 
rial. This simple model of the metamaterial interface has 
provided new insights into the non-Maxwellian boundary 
conditions that exist at the boundaries of metamaterials 
with spatial dispersion. 
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Appendix A: Calculation of s-polarized constitutive 
parameters 

Ref.© presents a simple and surprisingly successful one 
dimensional metamaterial homogenization model based 
on the same principle of replacing layers of a metamate- 
rial crystal with metasurfaces that was used in this paper. 
However, the formulas presented in Ref. |9| are for a p- 
polarized plane wave (E = E^e^ and H = H 2 e z ), while in 
this paper we consider s-polarized plane waves (E = E z e 2 
and H = Hj,ej,). The calculation of the constitutive pa- 
rameters for an s-polarized wave are very similar. Most 
of the changes can be made with the electromagnetic du- 



ality transformation F, z 



H», H„ 



fJ-b, 



and /ib —> though care must be taken with the new 
definitions of the field amplitude coefficients. We now 
briefly derive equations for the s-polarized one dimen- 
sional metamaterial homogenization model. 

The one dimensional homogenization model is a one 
dimensional system with all electromagnetic fields har- 
monic in time with frequency uj. The one dimensional pe- 
riodic spatial domain is separated into unit cells of length 
a. In the center of each cell is a metasurface with surface 
polarizability a defined in Eq. (IT5|) from the scattering 
matrix elements of a single layer of the crystal to be ho- 
mogenized. The space between each metasurface is filled 
with a homogeneous material with a background (sub- 
strate) permittivity £(,, permeability /if, and impedance 
z b = \/ Hb/tb- For one dimensional s-polarized plane 
waves, the Maxwellian equations for the microscopic elec- 
tromagnetic fields in between the metasurfaces are 
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de z 
dx 

dhy 
dx 



uj I 

1 — Hbhy = ~ 

c < 



(Al) 



u> 

i-e&e 2 

c 



The metasurfaces interact with each other through free 
plane waves with frequency uj and wavenumber fc = 
y/e b fi b uj/c. Across each metasurface, the e 2 and h y fields 
are discontinuous according to Eq. (1131) . The surface po- 
larization of the metasurfaces which defines the disconti- 
nuity is equal to the surface polarizability times the local 
field strength, as described in Sec. |TTJ Finally, the en- 
tire system is driven with external electric and magnetic 
current 



The field strength on the right hand side of Eq. (fT3|) is 
the sum of two solutions, the inhomogeneous solution to 
Eq. (|Al|) that does not obey the correct boundary condi- 
tions and the homogeneous solutions to Eq. (IA1|) which 
corresponds to the free waves mediating the interaction 
between adjacent metasurfaces. We choose our coordi- 
nates so that x = at the location of the nth metasur- 
face. The local value of the inhomogeneous solution at 
the location of the nth metasurface is 



iK- 1 




where the matrix K is defined as 



(A6) 
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T \(u)t- 



(A2) 



where the current strengths J3 and I2 are chosen arbitrar- 
ily. The solution to Eq. (|A1|) with the external current 
in Eq. (|A2|) is 



IhJ e b [i b uj 2 /c 2 - k 2 . 
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(A3) 

The field coefficients a b and b n are the electric field am- 
plitudes of free plane waves propagating in the posi- 
tive and negative x-directions respectively. These plane 
waves mediate interactions between adjacent metasur- 
faces. The subscript n specifies the metasurface that 
the plane waves interact with. For example the n and 
n + 1 metasurfaces interact with the right and left mov- 
ing plane waves with amplitudes a n and b n respectively. 
These field coefficients are related to each other by the 
Bloch phase condition 



K = 



(e b uj/c k x \ 
I k x fi b uj/cl 



(A7) 



The local value of the homogeneous solution at the 
location of the nth metasurface is defined as the average 
value from both sides of the sheet and is related to the 
field coefficients by 
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Combining Eq. (US]) with Eqs. (|A5IIA8|) . we find an equa- 
tion relating the field coefficients a n and b n to the exter- 
nal current. 
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which can be solved for the field coefficients to give 



, (A9) 






(A4) 



The field coefficients are also determined by the bound- 
ary conditions across each metasurface given by Eq. (|13p . 
The left hand side of Eq. (fT3|) is related to the field coef- 
ficients by 
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As in Ref.HI, the macroscopic polarization density is equal 
to the polarization of a metasurface divided by the unit 
cell length a. Using Eq. (|13l) . we can relate the field 
coefficients to the macroscopic polarization yielding 
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Here we have identified the matrix Q, defined in Ref. ;9; 
as the matrix relating the polarization and the external 
current. Incidentally, setting the determinant of —a + 

\A ■ B^ 1 equal to zero is an alternate way of deriving 
the dispersion relation in Eq. (fT9|) . Using Eqs. (|A5|) and 
(|A8|) . the matrix AB -1 is evaluated to be 

-2i /sm(k Q a)/z b sin(k x a) \ 

cos(k x a) + cos(fc a) I sin(A; x a) z b sin(k a)J ' 

(A12) 

Putting this into Eq. (|A11|) allows us to evaluate Q which 
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we can then use to evaluate the bulk susceptibility x 
defined in Ref. as 

X(u}, k x ) = (l - cj/cQ ■ K- 1 ) ■ Q. (A13) 
X is related to the constitutive tensor by 

C(u,k x ) = h M = t °) +x(w,k x ). (A14) 
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